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Abstract 

The Iterated Perturbation Theory (IPT) equations of the Dynamical Mean Field 
Theory (DMFT) for the half-filled Hubbard model, are solved on nearly real fre- 
quencies at various values of the Hubbard parameters, U, to investigate the nature 
of metal-insulator transition (MIT) at finite temperatures. This method avoids the 
instabilities associated with the infamous Pade analytic continuation and reveals 
fine structures across the MIT at finite temperatures, which can not be captured by 
conventional methods for solving DMFT-IPT equations on Matsubara frequencies. 
Our method suggests that at finite temperatures, there is an abrupt decrease in the 
height of the quasiparticle (Kondo) peak at a critical value of U c , to a non-zero, but 
small bump, which gradually suppresses as one moves deeper into the bad insulator 
regime. In contrast to Vollhardt & coworkers [J. Phys. Soc. Jpn. 74, 136, 2005] down 
to T = 0.01 of the half-bandwidth we find no T* separating bad insulator from a 
true Mott insulator. 
PACS: 71.30.+h 
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1 Introduction 



Feynman introduced a graphical representation of various perturbation schemes 
in the field theory which bears his namepQ, Feynman diagrams. In condensed 
matter physics, the natural high momentum cut-off provided by the lattice di- 
mensions prevents ultraviolet divergences. But yet remains the practical task 
of calculating them and interpreting the results. 
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Matsubara on the other hand introduced a very clever mathematical trick, 
which not only combines the thermal and quantum mechanical averaging re- 
quired in calculations of the correlation functions in a neat way and treats 
them on the same footing, but also provides a very convenient method for 
numerical calculations of the diagrams. The essential ingredient is that the 
Wick rotation r — it replaces the oscillatory e~^ kt factors by decaying e~^ kT 
factors in imaginary time, which are very convenient for putting on computer 
and convergence is particularly fast in iterative or self-consisted formulations 
of perturbation theory. However, the price will be payed when one has to undo 
the Wick rotation at the end of calculation to obtain dynamical quantities by 
replacing iu n u + irj, where r\ is an infinitesimal positive constant. 

The hurdle one faces in undoing the Wick rotation is that, if using the Pade 
approximation one fits a quotient of two polynomials /n(z) and quiz) to 
the table of data obtained for Matsubara frequencies iu n , and then replaces 
z — > uj + %r\ in the resulting function, the calculated spectral weight are not 
always stable with respect to variations in parameters N, M. Even if for some 
parameter regime, or for some particular problem one obtains relatively stable 
results, the aforementioned disrepute of the Pade approximation warns us 
about the reliability and/or the quality of the dynamical quantities obtained 
in this way. 

People have been worried about this question from time to time, and there 
has been some proposals for reliable way of using the Pade approximation 
in analytic continuation of numerical data: Beach and coworkers [2] proposed 
a symbolic computer aided algebra with arbitrary precision (typically 100- 
200 decimal places, which lack in single or double precision arithmetics of 
standard programming languages like C++ or Fortran). They also proposed 
a qualitative measure of the reliability of continued data. Mishchenko and 
collaborators [5fl] proposed an stochastic optimization method which allows 
one to handle both broad and sharp features of the spectrum on equal footing. 

On the other hands, Schmalian and colleagues proposed an alternative method 
which is quite intuitive and general [5]: Instead of solving the diagrammatic 
equations for Matsubara frequencies iu n , solve them for frequencies u + ij, 
where 7 is a finite constant. The finite value of 7 (usually taken to be less 
than the first nonzero Matsubara frequency) is important, and provides the 
damping required for convergence of the iterative solutions. At the same time 
analytic continuation from u + 27 to uj + ir), where 77 = + , not only is sta- 
ble, but also sustains fine features of the spectral function, such as shadow 
bands of the high temperature superconductors [5], as they found by applying 
this method to solve the diagrammatic equations of the fluctuation exchange 
approximation [7] . 

Dynamical mean field theory approximation has been successful in addressing 
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the issue of metal-insulator transitions in correlated electron systems at zero 
temperature [8]. Keeping the aforementioned concerns about the reliability of 
the Pade analytic continuation procedure in mind, in this paper we use the 
method of Schmalian et. a/[5] to re-examine the nature of MIT in the half-filled 
Hubbard model at finite temperatures. We find that paying the price at the 
beginning and solving slightly more difficult equations for u + 27 pays off and 
in addition to providing reliable Pade analytic continuation with absolutely 
no negative spectral weights, reveals fine structure in the insulating side of 
the MIT. A small bump in the spectral weight which persists in the insulating 
phase, is stably produced in our approach and can not be captured by Pade 
analytical continuation of the solutions of DMFT equations for Matsubara 
frequencies. 

The only difficulty with this method is that with conventional double precision 
numeric accuracy at lower temperatures it is hard to achieve convergence. But 
the asymptotic behaviors at T — > limit in our approach agree with other 
methods of solving the DMFT equations. 

The paper is organized as follows: First we analytically continue the IPT 
equations of DMFT to uj + 27 line above the real frequency axis. Then we 
present the numerical solutions of the resulting equations for various values of 
the Hubbard parameters, U, at half-filling, and elevated temperatures. Finally 
we present our conclusions. 



2 Formulation 



Within DMFT approximation, the problem of interacting electrons on a lat- 
tice can be mapped onto an effective impurity problem surrounded by a self- 
consistent bath. The impurity Green's function, Q, is related to it's bare coun- 
terpart via the Dyson equation [H], 



Go 1 = S + 

D(iu n + fi — E) 



where, 



D(iu n + v-X) = I de-. — (2) 

J %uj n + u — h — e 
—00 

is the Hilbert transform of density of states (DOS). In ([2]), D{iuo n + \i — E) 
is the on-site full Green's function for site o i.e. G OQ and it's imaginary part 
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;ives the interacting DOS, 

G(iu n ) = D(iu n + T*(iuj n )) (3) 



In IPT approximation, the self-energy is given by the second order perturba- 
tion theory [8] as, 



E(icu n )~U 2 J dr e^o(r) 3 . 
o 



(4) 



E( 2 ) 



DMFT equations written in Matsubara form yield no dynamical quantities, 
until the analytical continuation to real frequency axis is done, 

iuj n — > uj + ir) 



To see where lies the root of numerical problems, one notes that the real- 
frequency and imaginary time Green's functions are connected by [9] 

oo 

1 r e~ TUJ 

G(r) = - J du i + e _ u/T ImG{u + *0+) (5) 



where, G(t) = T J2n e JtJnT G(iuo n ) is the Fourier transform of Matsubara func- 
tion. Due to exponential factors, the small changes in G(t) (equivalently in 
G{iuj n )) are associated with large changes in G(u + irf). 

Now we turn our attention to the question of analytic continuation of DMFT 
equations in IPT approximation, parallel to the work of Schmalian and coworkers [H]- 
We rewrite the equation to be solved for nearly real frequencies u + ij, with 
a finite 7. Then, we go to the limit 7 — ► + via Pade approximation. Pade 
approximation at this stage turns out to be stable. The finite parameter 7 is 
chosen to provide the attenuation factors (as will be seen below) needed for 
convergence of the self-consistent equations. 

The first equation to be continued analytically to nearly real axis is (j3J) , which 
bears no difficulty, 

G{u + i-f) = D(u + i-f + n-Z(Lj + i-f)). (6) 
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The next equation to be continued to nearly real-frequency axis is ([T]), which 
again simply reads 

Go V + il) = E(w + il) + - 1 (7) 

L>(cj + ry + /i - L(cj + 27)) 



The main problem relies on analytically continuing the IPT approximation, 
Eq. to nearly real frequency axis. This equation depends on the frequency 
not only through the Fourier factor e luJnT , but also through the Green's func- 
tion &o(t). In this case, analytical continuation must be performed through 
the change of integral to it's retarded form. Here the difficulty arises from the 
fact that we've solved the action equation (for a detailed review see Georges et 
a/.[H], section III. A), on the imaginary axis to yield Matsubara function, so if 
one wishes to gain the physical quantities (such as retarded green's function) 
one must either solve the problem originally on real frequency axis or analyt- 
ically continue it. Kajueter and Kotliar [TOlTT] made an ansatz for self-energy 
on the real frequency axis of the form 



where, E^ 2 -* is the second order contribution to self-energy from (j3J). Of course, 
this equation alone doesn't solve the problem, since there are other functions 
written in Matsubara form. Fortunately there is another way to overcome the 
problem. Equation (j4]) can be written as 

P 

Z{iu n ) ~U 2 J dr e l ^ T go(r)g (r)g (-r). (9) 




Using the Fourier transformation 

00 

0o(t)= £ e^ T g (iu n ), (10) 

n=— 00 

it can be written as 

P 

fdTe^ T g (r)g (r)g (-T)=Y / G ^ k )x ^(u; n + u k ))^ (11) 







k 



where, 
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X°(i(Wn + = J] <7o(«^)C?o(«(^n + Wj + Wfc)) 



(12) 



is the particle-hole bubble. 



With the aid of contour integration and employing the complex forms of Fermi 
and Bose functions, 



/(*) 



1 



n(z) 



1 



3 f3z _ X 



(13) 



which have their poles exactly at Matsubara frequencies z = iu n = i(2n + 
l)7r//3 and z = iu' n = i(2n)n/j3, respectively, the summation over imaginary 
frequencies can be done, 



i 

de 
2T\i 



(14) 



{f(e + %i)g G {wk + e + ii)Qo{e + %i) 

- f(s - ii)Go{iv k + £ - il')Go(£ + h') 
+ f(e + ii)Q Q {e + ii)Q Q {e + ij - iv k ) 

- /(e - <Y)</o(e - iV)0o(e - «Y - «ffc)}- 

Now, performing the analytical continuation iz>% — * u; + ry, and making use of 
Kramers-Kronig transformation for the green function, along with the Laplace 
transformation 



1 



(15) 



the particle-hole bubble reduces to 



o 



(16) 



where, 



X (t) = -z(2tt) 2 (p(t) [A*(t)+A(-t) 



p*(t) + A*(-t)e 27 * 



(17) 
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In the above formula the spectral densities are given by, 
1 1 



P(t) = -o / d£ Im ^(e + il)e~ l£ \ (18) 



and, 



oo 

I 



A{t) = — / def*(e + ii)QZ(e + t 7 )e— (19) 

Z7T 



-OC 



In an analogous manner, the Laplace transform of self-energy can be written 
as: 



£(t) = -ilixT U 2 Re X °(0 + i0 + ) e*p(t) 

+i(2ir) 2 U 2 x (v°(t)[A(t) + A*(-t)e 27 *] 



-p(t)[5*(t) + 5(-t)e^]) (20) 

where, 



1 1 



A*) = -7T- I d£ Im X°( £ + ^) e ~ l£t > ( 21 ) 

Z7T 71 J 



and 



oo 

B{t) = — ( den*(e + i-f) X °*(e + i-f)e- l£t . 



[22] 



3 Numerical Results 



The numerical results are given in units of the half bandwidth W/2 = 1. The 
temperature T = 0.048 and the Hubbard parameter U varies from U = 
through 7.5. 

To perform the Fourier transforms we used the FFT routines. To overcome 
the alias effects arising from the 1/uo tail, we chose the energy mesh large 
enough to ensure that the function vanishes at the boundaries. We have used 
N = 2 15 frequency points in an energy range of 2u max = 250. Also we have 
taken 7 = 7rT/4. The finite value of 7 provides necessary damping factors in 
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the convergence of the iteration process. Therefore at lower temperatures it 
becomes harder to achieve convergence. Convergence at lower temperatures re- 
quires smaller values of 7 and hence higher accuracy than conventional double 
precision [2]. 

There are many lattices examined by the DMFT method, some of which is 
mentioned in the review paper by Georges et al. [S|. We work with Bethe 
lattice with semicircular DOS [T2j . 

We begin with an initial guess for the green function, Qo{uj + £7). Using 
Eqs. f fT8][T9l) we calculate p{t) and A(t) which are used to evaluate X°(t) in 
Eq. ([17]). Then we calculate v°(t) and B(t), Eqs. (ETHTSi . through which the 
self-energy £(£), Eq. (I2U|) . and hence its Laplace transform, E(u; + ry), can 
be evaluated. To close the iteration loop, E(u; + £7) is inserted in the self- 
consistency equation ([1]) to update the initial Greens' function and iteration 
is performed until convergence is achieved. 

We have obtained converged solutions for U in the range [0, 7.5] for T = 
0.048. At this temperature the transition occurs at critical value U c ~ 3.0 + , 
in agreement with the results of Zhang and coworkers [13], and Vollhardt et 
al. [H. 

3.1 The spectral density before the transition 
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Fig. 1. DOS at T = 0.048 and values of U < U c in the metallic side. When the 
Hubbard parameter is increased the weight of fine peak is transferred to it's neigh- 
bors, unlike the zero temperature case, the height of the central Kondo peak is not 
constant. 

Before the transition in the metallic regime, as one increases U the DOS tends 
to split, leaving a central Kondo peak corresponding to quasiparticle at the 
Fermi level (Fig. [1]) in agreement with the conventional methods of solving 
the DMFT equations [8] • In our computer code we have not forced the height 
of the Kondo peak to be constant. In our finite T solution unlike the zero 
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temperature case the hight of the Kondo resonance is not constant during 
the transition [15]. The height of the quasiparticle peak at Fermi energy is 
instead gradually redistributed and shifted to the upper(lower) edge of the 
lower (upper) Hubbard band. 

3.2 The spectral density near the transition 

The transition occurs at a critical Hubbard parameter U c . We find the value 
for this critical Hubbard parameter in the range 2.99 < U c < 3.01. At the 
transition the height of peak falls off suddenly to a very small, but nonzero 
value, Fig. [2j 

1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 


-2-1012 
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Fig. 2. The evolution of DOS when the transition occurs. The results is shown for 
two representative values U = 2.9 and U = 3.1 at temperature T = 0.048. 



3.3 The spectral density after the transition 

The evolution of the spectral function after the transition is shown in Fig. [3] 
In our solution increasing U results in depletion of the spectral weight but 
does not give a clean gap at finite temperatures. For larger values of U the 
spectral density at Fermi level is still finite and vanishes only in the limit 
U — > oo, so that there is no real transition but a crossover from a metallic-like 
to an insulating-like solution. Vollhardt et a/.[H] named this state as a bad 
insulator. 

The finite temperature analogue of the Mott-Hubbard gap can be thought of 
as the peak-to-peak distance between the upper and lower Hubbard bands 
in Fig. [31 As can be seen this distance increases proportional to U. With 
increasing U the "gap" tends to expand and the spectral weight in this region 
tends to decrease, but as we mentioned the gap is not clean at T = 0.048. At 



U=2.9 
U=3.1 
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Fig. 3. DOS after the transition, at temperature T = 0.048. Increasing U tends 
to push upper and lower Hubbard bands far apart, but the little peak resits and 
remains there at even much larger values of U. 

U — > oo one expects to obtain a clean gap around the Fermi level. In order to 
obtain clean gap at finite large U one has to go to T — > limit. 

In the solution obtained by Vollhardt et al. [33] they introduce a temperature 
T* below which the bad insulator becomes a true insulator. Down to T = 0.01 
where we have obtained converged solution in our method we do not find such 
a T*. The trend of the lower T data in Fig. H] shows that at T — > limit 
the gap starts to become cleaner and develops sharp edge characteristic of 
IPT method [S|. The small central peak in the bad insulating side gradually 
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Fig. 4. The evolution of the spectral weights in the bad insulator regime for U = 5.0 
for lower temperatures. The lower the temperature the cleaner the gap and the 
insulator evolves to Mott insulator. As the figure shows, the walls of the gap tends 
to vertical direction at T — > limit. 

suppresses with decreasing T but does not vanish at any finite U (Fig. [4]). 
The filling of the Mott-Hubbard gap by spectral weight transfer with increas- 
ing temperature has recently been detected experimentally by photoemission 
experiments [IE] . 

An important aspect of insulating phase obtained in the present method is 
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the presence of a small bump at the Fermi level which is clearly seen in the 
Figs. [3] and HI This little peak has a height about 10% of neighboring bands 
and does not vanish after the transition at any U. 



4 Summary and conclusions 

We employed a new method for solving diagrammatic equations to solve IPT 
equations of DMFT. We first analytically continued the equations to nearly 
real frequencies. Solving the continued equations is slightly harder, but seems 
to capture fine features of the spectral weight not previously reported. 

We found that the height of Kondo peak does not remain constant before tran- 
sition which is in contrast to the zero temperature result of Mixller-Hartmann [T3] , 
which states that the height of the Kondo peak must remain constant before 
the transition occurs. The transition is first order and the Kondo peak sud- 
denly falls off to a very small but non-zero value. For larger values of U this 
small spectral density at Fermi level will persist and vanishes only in the limit 
U — > oo or T — > 0. In this limit the gap also becomes cleaner and we obtain 
a true insulator only in T — > limit. Hence we do not find a T* [T41. 



Acknowledgement 

M.B.F. would likes to thank Prof. M.A. Vesaghi for his many suggestions and 
constant support during this research. S.A.J, was supported by ALAVI Group 
Ltd. 



References 

[1] H. Bruus, and K. Flensberg: Many-body Quantum Theory in Condensed Matter 
Physics, Oxford University Press, (2004). 

[2] K.S.D. Beach, R.J. Gooding and F. Marsiglio: Phys. Rev. B 61 (2000) 5147. 

[3] A.S.Mishchenko, N.V.Prokofev, A. Sakamoto, and B.V. Svistunov: Phys. Rev. B 
62 (2000) 6317. 

[4] A.S. Mishchenko, N. Nagaosa, N.V. Prokofev, A. Sakamoto, and B.V. Svistunov: 
Phys. Rev. B 66 (2002) 020301. 

[5] J. Schmalian, M. Langer, S. Grabowski, K. H. Bennemann: Computer Physics 
Communications 93 (1996) 141. 



11 



[6] M. Langer, J. Schmalian, S. Garbowski, K.H. Bennemann: Phys. Rev. Lett. 75 
(1995) 4508. 

[7] N.E. Bickers, D.J. Scalapino: Ann. Phys. (NY) 193 (1989) 206. 

[8] Georges and G. Kotliar and W. Krauth and M. J. Rozenberg: Rev. Mod. Phys. 
68 (1996) 13. 

[9] G. Dopf et al: Helv. Phys. Acta 65 (1992) 257. 

[10] H. Kajueter and G. Kotliar: Int. J. Mod. Phys. 11 (1997) 729. 



[11] K. Held: |arXiv:cond-mat/0511293| (2005). 

[12] E. N. Economou: Green's Function in Quantum Physics (1983) 2nd Ed. 

[13] X.Y. Zhang, M.J. Rozenberg, and G. Kotliar: Phys. Rev. Lett, 70 (1993) 1666. 

[14] D. Vollhardt, K. Held, G. Keller, R. Bulla, Th. Pruschke, I. A. Nekrasov and 
V. I. Anisimov: J. Phys. Soc. Jpn. 74, 136, 2005 

[15] E. Muller-Hartmann: Int. J. Mod. Phys. B 3 (1989) 2169. 

[16] S. K. Mo, H. D. Kim, J. W. Allen, G.-H. Gweon, J. D. Denlinger, J. H. Park, 
A. Sekiyama, A. Yamasaki, S. Suga, P. Metcalf, and K. Held: Phys. Rev. Let. 
93 (2004) 076404. 



12 



